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We report the results of molecular dynamics (MD) simulations of a complete monolayer of 
hexane physisorbed onto the basal plane of graphite. At low temperatures the system forms a 
herringbone solid. With increasing temperature, a solid to nematic liquid crystal transition 
takes place at T\ = 138 ± 2K followed by another transition at T 2 = 176 ± 3K into an 
isotropic fluid. We characterize the different phases by calculating various order parameters, 
coordinate distributions, energetics, spreading pressure and correlation functions, most of 
which are in reasonable agreement with available experimental evidence. In addition, we 
perform simulations where the Lennard- Jones interaction strength, corrugation potential 
strength and dihedral rigidity are varied in order to better characterize the nature of the 
two transitions through. We find that both phase transitions are facilitated by a "footprint 
reduction" of the molecules via tilting, and to a lesser degree via creation of gauche defects 
in the molecules. 

PACS numbers: 64.70.-p, 68.35. Rh, 68.43.-h 



I. INTRODUCTION 

The physical adsorption (physisorption) of 
atoms and molecules onto a substrate to form 
a quasi two-dimensional (2D) film has resulted 
in the observation of a rich variety of behavior 
not realized in the system's corresponding three- 
dimensional (3D) or "bulk" state. 

Among the species used as adsorbates, alka- 
nes are of considerable interest because of their 
technological importance (e.g. as lubricants) 
and because they are among the simplest fami- 
lies of molecules of compounds whose members 
differ mainly in their length. In fact, straight- 
chained n-alkanes represent a fine balance be- 
tween complexity (i.e. many internal degrees of 
freedom which also affect the interaction between 
molecules) and the simplicity of their structure 
(as compared to other organic molecules). In 
fact, over the past decades there has been a 
constant and renewed effort to better under- 
stand the adsorption of alkaknes over a variety 
of substrates. Among the substrates used for 
physisorption studies, graphite, with its excel- 
lent mechanical stability, availability and sym- 
metry, has proven to be one of the best choices. 
A considerable amount of experimental, theoret- 



ical and computational work has, thus, been de- 
voted to these studies. 

Hexane on graphite was first studied by Krim 
et al. using LEED as well as neutron diffraction. 
For submonolayer coverage a uniaxial incom- 
mensurate (UI) herringbone phase is observed 
below ca. T = 151 K, where a first-order melting 
transition is found (but note that but a determi- 
nation of the molecular orientations was not pos- 
sible). At low temperatures, as the coverage is 
increased, the UI phase evolves continuously into 
a 2x4\/3 commensurate structure at completion. 
More recently, Taub and co-workers 2 completed 
extensive neutron 2 ' 3 and X-ray 2 diffraction stud- 
ies of hexane on graphite for submonolayer 2 , 
monolayer 2-6 and multilater 2 coverages. Their 
findings indicate that at low temperatures a 
complete monolayer forms a commensurate her- 
ringbone structure which evolves with increas- 
ing temperature into a rectangular centered 
solid/liquid coexistence region by ca. 150 K and 
melts at around 175 K 2 ~ 6 . For submonolayer 
coverages, these authors proposed a structure 
corresponding to a UI phase comprised of com- 
mensurate regions separated by low-density fluid 
filled domain walls 2 . 
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In the 1990's, Hansen et al. reported stud- 
ies combining MD computer simulations with 
neutron and x-ray diffraction studies of bu- 
tane and hexane on graphite at monolayer 
completion 4,5 . Experimentally, butane melts 
abruptly at around T = 116 K directly from 
a solid commensurate rectangular-centered her- 
ringbone (HB) phase into a liquid. Orientational 
ordering about the surface normal is lost through 
rotation about the center of mass, which is si- 
multaneous with melting. Hexane, on the other 
hand, undergoes a loss of translational order at 
about 150 K into a phase with short-range order 
which is thought to involve mobile rectangular- 
centered (RC) islands within a fluid. Then, at 
around T = 175 K the system melts into a fluid. 

The computer simulations of Hansen et aZ. 4 ' 5 
showed that, for hexane, the internal formation 
of gauche defects coupled with out-of-plane tilt- 
ing is concurrent with melting. The necessary 
creation of in-plane room by either gauche de- 
fects or tilting at the onset of melting is referred 
to as a "footprint reduction" mechanism. These 
simulations showed, furthermore, that the na- 
ture and temperature of the transition is very 
sensitive to the presence of gauche defects, as 
the melting temperature raises drastically when 
gauche defect formation is suppressed. However, 
the resulting melting temperature of a hexane 
monolayer was found to be 222 K, significantly 
different from the experimental value of 175 K. 
This discrepancy was addressed in a series of MD 
simulations of hexane on graphite by Velasco 
and Peters 7 which found a better agreement with 
experimental values of the melting temperature 
when a considerably lower adsorbate-interaction 
strength was used, and under those conditions 
molecular gauche defect formation becomes al- 
most irrelevant to melting. (In this paper we find 
the presence of an intermediate nematic phase. 
In the herringbone solid to nematic phase tran- 
sition we find gauche defects to be irrelevant. 
For the nematic to isotropic liquid phase tran- 
sition, we find that both mechanisms contribute 
to the "footprint reduction", although it is ap- 
parent that tilting is the dominant effect.) 

More recently, Peters and coworkers have 
also examined the behavior of monolayer 8 and 
bilayer 9 hexane on graphite using MD simulation 



techniques. The monolayer study compared MD 
results based on the standard (isotropic) unified 
atom (UA) approximation to simulations where 
an anisotropic force field was employed. Using 
a computational cell of 112 hexane molecules, 
and performing MD simulations corresponding 
to 200,000 steps of 1 fs both models yielded 
the same basic physics, especially in regards to 
the characterization of the phases observed (the 
transition temperature from commensurate her- 
ringbone to the rectangular centered orientation- 
ally ordered phase is found at around 150 K and 
the melting temperature at 175 K in both cases, 
the main difference between the two models is 
that the anisotropic model promotes more in- 
plane mobility from vacancy creation due to in- 
creased molecule tilting at all temperatures). 

Even though hexane on graphite has been 
studied quite well experimentally, it has not 
been as thoroughly investigated computation- 
ally, and there is still some doubt about the 
roles played in the phase transitions by the inter- 
nal degrees of freedom of the molecules (mainly 
gauche defects) and their relative importance as 
the adatom-substrate interaction is varied. Fur- 
thermore, it is not self-evident that the simu- 
lations have achieved ergodicity to a satisfac- 
tory degree, as only a relatively limited simula- 
tion time has been explored so far (maximum 
of 2 ns in reference 8). The general purpose 
of the work reported here is to enrich our un- 
derstanding of phase transitions in physisorbed 
molecular systems, in particular of hexane in 
graphite. Specifically, we wish to: (i) more 
sharply delineate the phase transitions of hex- 
ane on graphite using longer simulations (up to 
5 times longer than previously reported); (ii) to 
understand the effects of varying the graphite 
corrugation strength and the adatom-adatom in- 
teraction strength on the observed phase tran- 
sitions; (Hi) to further delineate the role of 
the internal degrees of freedom of the adsorbed 
molecules, i.e. the bond bending and dihedral 
torsion in the observed phase transitions, and 
(iv) to gain insight in the energetics and possible 
mechanisms for the observed phase transitions. 

Since the system under consideration is, to 
some degree, (quasi-) two-dimensional, one may 
have expected that the observed phase transi- 
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tions would show some of the universal features 
expected from the beautiful theory of 2D melt- 
ing developed in the 1970's by Kosterlitz, Thou- 
less, Halperin, Nelson and Young 10-14 . However, 
the importance of the out-of-plane behavior of 
the hexane monolayer, and the "polarizing" ef- 
fect of the substrate corrugation make it diffi- 
cult to reliably place these transitions in a gen- 
eral theoretical context. Furthermore, tt is ev- 
ident that the the relatively small system sizes 
achievable in MD simulations does not permit 
a complete characterization of the phase tran- 
sitions observed (e.g. their order), but some- 
times the modality of energy distributions near 
the critical regions may give indications on what 
behavior to expect. 

The paper is organized as follows. In Sec. II 
we discuss the computer simulation procedures. 
The results of the simulations are presented in 
Sec. Ill, and are discussed in detail in Sec. IV. 

II. COMPUTATIONAL DETAILS 

The hexane molecules in our study are mod- 
eled with two methyl (CH3) pseudo atoms on its 
ends and four methylene (CH2) pseudo atoms 
in between. The lowest energy conformation of 
the molecule consists of all pseudo atoms being 
co-planar, forming a ziz-zag pattern with 114° 
three-body bond angles (dihedral angles) and 
bond lengths (nearest neighbor separations) of 
1.54 A. In the simulations that follow we keep 
the the bond lengths fixed, but allow the other 
internal degrees of freedom to vary (bond an- 
gle bending and dihedral angle torsion). The 
graphite substrate modeled as being static, con- 
taining an infinite number of stacked graphene 
sheets (see Fig. 1). The structure contains point 
carbon atoms at the vertices of all bonds and the 
bond length is taken to be 2.46 A. 

A. Method 

A constant particle number, planar density 
and temperature (N = 672, p = 1, T) molecu- 
lar dynamics (MD) method is employed to con- 
duct atomistic simulations a system of N rn = 
2 x m x n = 112 hexane molecules at monolayer 
completion (m = 4 and n = 14). The rectan- 



gular computational cell (see Figure 1) has its 
corners coincident with graphite hexagon cen- 
ters and has dimensions m x 4V3a ~ 68.1735 
A and n x 2a = 68.88 A in the x and y di- 
rections respectively. Such a cell gives a planar 
density of p = 0.02385 molecules/ A 2 and an area 
per molecule of 42 A 2 . Each hexane molecule 
consists of six united atom (UA) pseudoatoms- 
two methyl (CH3) and four methylene (CH2). 
Periodic boundary conditions are utilized for 
each pseudoatom in the x, y directions and free 
boundary conditions are present in the vertical 
z direction. To perform the simulations at con- 
stant temperature, we periodically rescale the 
velocities so that equipartition is satisfied for the 
center-of-mass, rotational and internal tempera- 
tures: 
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respectively (no significant differences were 
found when "thermalizing" a subset of the tem- 
peratures above). Here N m is the number of 
molecules, nc the number of pseudoatoms per 
molecule, the i index runs over the molecules 
(1 < i < N m ) and the j index runs over the pseu- 
doatoms within a given molecule (1 < j < nc)- 
Furthermore, for each molecule, M, is the mass, 
Vi,cM is the center-of-mass velocity, ui is the 
angular velocity and 2j is the moment of iner- 
tia tensor (relative to the COM). Last, is 
the velocity of the j-th pseudoatom in the i-th 
molecule, and k B is Boltzmann's constant. 

The time step for all simulations is cho- 
sen to be 1 fs and integration of the equa- 
tions of motions is achieved with a velocity Ver- 
let RATTLE 15 algorithm which keeps the pseu- 
doatom bond lengths constant. Runs are typi- 
cally started from a herringbone configuration, 
as suggested by experimental evidence (e.g. the 
neutron diffraction data of Ref. 1), and allowed 
to equilibrate for roughly 3-5 xlO 5 steps, aver- 
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T=130K 



T= 160K 



T=200K 



FIG. 1: Simulation cell showing snapshots of typical (equilibrium) configurations for the hexane/graphitc 
system in the commensurate herringbone solid (T = 130 K), the orientationally ordered nematic (T = 160 
K) and isotropic liquid (T = 200 K) phases. CH 3 end-of-chain pseudoatoms are red, CH 2 are yellow 
and the C-C bonds in the graphite substrate are purple (the graphite lattice spacing is a = 2.46 A). The 
computational cell has dimensions m x 4V3a ~ 68.1735 A and n x 2a = 68.88 A in the x and y directions 
respectively, where m = 4 and n = 14, and periodic boundary conditions are used in x and y directions. 
There are 2 x m x n = 112 hexanc molecules. 



ages are taken for 0.5-1.0xl0 6 steps after equi- 
libration. Some runs are started from the final 
configuration of some lower temperature run. A 
good measure of convergence is that incorporat- 
ing the results from both types of initial condi- 
tions does not create noticeable scatter in the 
results. The only hysteresis observed is that, as 
is usual in these type of simulations, "freezing" 
by starting from a high temperature phase and 
lowering the temperature is not readily achieved. 

B. Interaction Potentials 

The model adopted for our study of adsorbed 
hexane molecules on graphite corresponds to 
having both non-bonded and bonded interac- 
tions. The first of two non-bonded interactions 
is the adsorbate-adsorbate interaction between 
pseudoatom atom i and j separated by a dis- 
tance rjj, which is modeled by a Lennard- Jones 
pair potential 
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are applied in order to describe mixed interac- 
tions when particles i and j are of different types; 
lattice sums are taken out to 10 A and values for 
the potential parameters used in our work 16 are 
given in Table I. 

The second non-bonded interaction is the 
adatom-graphite interaction. The graphite is 
modeled as being of infinite extent in the xy 
plane and semi- infinite in the vertical — z direc- 
tion. The adatom-graphite interaction potential, 
however, is not taken as a discrete sum over the 
carbon atoms. Steele's expansion 17 is utilized 
because it exploits the fact that the symmetry 
of the substrate breaks the adatom-graphite in- 
teraction into a vertically (laterally averaged) 
portion and a much weaker portion related to 
the undulation in the interaction due to the 
graphite hexagonal symmetry. Such a represen- 
tation saves a considerable amount of computing 
time. In other words, the interaction of the ad- 
sorbate atoms with the substrate is the Fourier 
expansion 
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TABLE I: Non-bonded potential parameters used in 
the simulations 



Parameter 


Value 


e CH 2 ; e CH 3 


72 K 


CTCH 2 , fCHg 


3.92 A 


q 


2 


a 


2.46 A 


a s 


5.24 A 2 


d 


3.357 A 




44.89 K 


(7 gr 


3.66 A 



Here g n is the modulus of the n-th graphite re- 
ciprocal lattice vector and the K n (x) are modi- 
fied Bessel functions of the second kind. The in- 
teraction is obtained by summing over an infinite 
number of graphene sheets. Only fi(xi, yi) is de- 
fined above because the sum in x, y-dependent 
part of equation (4) converges extremely rapidly 
and only the n = 1 term is needed. All potential 
parameters for the adsorbate-substrate interac- 
tion are identified in Table I. 

There are two types of bonded interactions 
considered in this work: bond angle bending 
and dihedral angle bending (we assume that the 
bond lengths are fixed at 1.54 A, the RATTLE 
algorithm allows for constrained solution of the 
equation of motion 15 ). The bond angles are as- 
sumed to be harmonic; their bending potential 
energy comes from Martin and Siepmann 16 and 
is given by 
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bond angle and kg is the angular stiffness. In ad- 
dition, the expression for the energy of dihedral 
bending (torsion) 8,18 is of the form 
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where 
stants. 



d is the dihedral angle and the Cj are con- 
All values for the bonded potential pa- 
rameters are shown in Table II. This potential 
has one global minimum for the trans configu- 
ration (fid = and two local minima with ~234 
K higher in energy for the two gauche configu- 
rations at (fid = ±2-7r/3. 

The full potential energy can thus be written: 



u = u nh + u h . 
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TABLE II: Bonded potential parameters used in the 
simulations 



Parameter 


Value 


kg 


62,793.6 K/rad 2 


Oo 


114° 


CO 


1,037.76 K 


Cl 


2,426.07 K 


C2 


81.64 K 


C3 


-3,129.46 K 


C4 


-163.28 K 


C5 


-252.73 K 



where Of, is the bond angle, Of, is the equilibrium 



III. RESULTS 

We performed MD simulations as described 
above for many temperatures in the 100-200 K 
range. Figure 1 shows snapshot of typical config- 
urations of the simulated hexane/graphite sys- 
tem for three characteristic temperatures where 
the herringbone (HB) commensurate solid, the 
"nematic," and the fluid phases are observed. 
It should be mentioned that promotion of 
molecules to higher layers was not very sig- 
nificant at any of the temperatures studied. 
From this figure it is evident that in the low- 
temperature phase molecules are arranged in a 
HB structure where the molecules' azimuthal 
angles fa £ [0,7r] 19 take predominantly values 
near {30°, 150°} {fa is defined as the angle that 
the axis of the smallest moment of inertia for 
molecule i makes with the x-axis). Moreover, 
the molecules are in registry with the substrate. 
Clearly this HB phase is solid. In the high- 
temperature regime, it is also quite clear that the 
molecules are randomly positioned and oriented 
(and in particular not in registry with the sub- 
strate), which is strongly suggestive of a liquid 
phase. The situation of the mesophase, however, 
is not as evident. It is clear that it has a defi- 
nite orientational order, quite different than the 
HB (most azimuthal angles centered around 0°, 
the molecules appear not to be commensurate 
with the substrate, and it is not clear whether 
the system is in a liquid crystalline state (i.e. a 
nematic), or in a solid phase (i.e. a rectangular 
centered solid). The exact nature of the different 
phases (and in particular of the mesophase) is 
addressed more quantitatively below, by means 
of examination of various order parameters, en- 
ergy analysis, spreading pressure, distributions 
and correlation functions. 

A. Order parameters 

Figure 2 shows the temperature dependence 
of several molecular order parameters which are 
useful in characterizing the system's behavior 
through and across phases. The herringbone or- 
der parameter is used to quantify the in-plane 
orientational ordering of the molecular axes and 
is defined as 
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FIG. 2: Structural order information. The follow- 
ing order parameters (OP) are presented as func- 
tions of temperature: Herringbone OPherr (solid cir- 
cles), center of mass OP C omm (open triangles), ne- 
matic OP nem (solid diamonds) and tilt OPuu (solid 
squares) . 



where the sum is peformed over all molecules, 
the integer j is assigned so as to take the two 
sublattices into account and the averages corre- 
spond to sampling a large number of uncorre- 
cted configurations in each MD simulation. For 
a "perfect herringbone" with fa £ {45°, 135°}, 
OPherr = 1] obviously OPherr vanishes if the two 
sublattices cannot be distinguished (e.g. if all 
angles are clustered around one angle or if all an- 
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gles are sampled equally) either statically (over 
the lattice) or dynamically (over time). For the 
HB phase, where the fa's are clustered around 
30° and 150° (for further evidence see the 130 
K curve in Fig. 7), we find OPherr — — 
0.866, as observed for temperatures smaller than 
Ti = 138 ± 2 K. Not surprisingly OP herr ~ for 
the high temperature phase as essentially all az- 
imuthal angles are equally sampled (see the 200 
K curve in Fig. 7). We also find OPherr — for 
the mesophase, which is also evident from the 
(pi distributions which show a single significant 
orientation and thus the lack of two angular dis- 
tinct sublattices (see the 160 K curve in Fig. 7). 

To investigate whether the adsorbed 
molecules are in registry with the graphite 
substrate we also measure the commensurate 
order parameter OP comm , defined as 



OP r 
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where the outer sum runs over all N m molecules 
at center-of-mass position rj, and the inner 
sum runs over all six graphite reciprocal lattice 
vectors g s £ {(±7r/3a, 0), (=F7r/6a, ±ir/y/3a), 
(±7r/6a, ±7r/ y/3a)}. OP comm takes on a value of 
unity were all molecular centers to remain stat- 
ically graphite hexagon centers and vanishes in 
the limiting case of the molecular centers uni- 
formly sampling positions in the xy plane. As 
with OPherr, Such a sampling could occur either 
in a static fashion (an infinite incommensurate 
solid) or in a dynamic fashion such as in a non- 
registered film. It is evident from Fig. 2, that as 
the HB phase disappears at T\ = 138 ± 2 K, the 
film also looses registration with the substrate. 
This is also evident in the decrease of the corru- 
gation energy (see Sec. IIIB and bottom panel 
of Fig. 3). 

Both OPherr and OPcomm are useful in pro- 
viding information on the solid to nematic tran- 
sition at T\ = 138±2 K but are unable to provide 
any characterization of the nematic to isotropic 
liquid transition. The nematic order parameter 
OPnem provides precisely the needed informa- 
tion and is defined by 
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where (pdir is an "average" azimuthal direction 
defined by maximizing the nematic order param- 
eter: dOPnem/d(f) d i r = 0: 
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Note that due to the 1/2 pre-factor, it is impor- 
tant to use the 4-quadrant version of the tan -1 
function to uniquely define <pd%r (the 2-quadrant 
version leads to the possibility of n/2 rotations 
which are relevant for the final value of OP ne m 
in some cases). 

The center panel of Fig. 2 shows the tem- 
perature dependence of the nematic order pa- 
rameter. This order parameter clearly indicates 
that, in addition to the HB to nematic phase 
transition at T\ = 138 ± 2 K, a second orienta- 
tional phase transition occurs at T2 = 176 ± 3 
K. As shown earlier for the OPherr, OP ne m also 
is clearly correlated to the azimuthal angle dis- 
tributions (Fig. 7). For a "perfect herringbone" 
solid (OPherr = l)j (ftdir is completely undefined 
and OPnem vanishes, but for the solid commen- 
surate herringbone phase of hexane on graphite 
(OPherr = v3/2), (frdir = and the limiting 
value of OP nem is 0.5. In a perfect nematic 
phase, all molecules share the same orientation 
and OP nem = 1- Our results of max[OP nem ] ~ 
0.8 show the presence of significant fluctuations 
in the molecular orientations and of the pres- 
ence of differently oriented domains, both fea- 
tures are also clearly qualitatively visible in the 
160 K data presented in the configuration snap- 
shot (Fig. 1) and azimuthal distribution (Fig. 7). 
In the liquid phase we would expect OP ne m to 
vanish, the small finite value in the high temper- 
ature phase is an indication of a small residual 
order likely due to the small system size. Note 
that contrary to the previous order parameters, 
OPnem does not vanish due to dynamic fluctua- 
tions leading to an overall rotation of the aver- 
age director over time. We would like to remark 
that our results for the nematic order parameter 
(center panel of Fig. 2) differ from those found 
by Peters and Tildesley (PT) 8 : (a) they find 
OPnem — —0.5 for the HB phase; (b) they find a 
relatively small max[OP nem ] ~ 0.0-0.4; (c) they 
find the HB to nematic transition at a higher 
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temperature (ca. 150 K); and (d) they find a rel- 
atively undefined nematic to liquid phase tran- 
sition at ca. 175 K. The difference in (a) (and 
to some extent (b)) is due to a factor of 1/2 
difference in the definition of </><& r in equation 
(12). The higher "importance" of the nematic 
mesophase in our work [points (b-d)] we belive 
is due to the considerably better statistics in our 
runs (whereas PT used 5 x 10 4 thermalization 
and 2 x 10 5 averaging steps, we typically ran our 
simulations for 3-5 xlO 5 thermalization and 0.5- 
l.OxlO 6 averaging steps). 

The tilt order parameter is designed to pro- 
vide a measure of the out-of-plane orientational 
order of the system— in particular, the amount of 
molecular tilting — and is defined as 



IPtat = — (^(P 2 (cos^ 
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2N n 
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were P2(x) is the Legendre polynomial and B{ 
is the angle that the smallest moment of iner- 
tia axis of molecule i makes with the normal to 
the substrate. Clearly OPtut is closely associ- 
ated with the tilt angle distribution P(B) (see 
Fig. 8). When all molecules are parallel to the 
surface OPtut = —0.5, and OPtut = 1 when all 
molecules are perpendicular to the surface. The 
bottom panel of Fig. 2 shows the temperature 
dependence of OPtut- At low temperatures the 
molecules clearly are (as expected) roughly par- 
allel to the graphite substrate to maximize their 
adsorption energy. At the HB to nematic tran- 
sition, OPtut has a significant peak, indicating a 
tendency of the molecules to librate our of the 
surface plane, but this tendency is reduced again 
as the nematic phase is better "established." 
For increasing temperatures OPtut barely in- 
crases until the nematic to liquid phase transi- 
tion, where it starts a significantly faster growth. 
Tilting of the molecules can obvioulsy be a signif- 
icant factor in the phase transitions, as the "foot- 
print" of the molecules is reduced. To put things 
in perspective, however, OPtut (200 K) ~ —0.44 
roughly corresponds to an average libration of 
only 11° out of the plane. Nevertheless, for a 
full monolayer coverage the resulting "footprint 
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FIG. 3: Thermodynamic information. Ensemble av- 
erages per molecule of the Lennard- Jones interac- 
tion energy (Ulj) and the lateral corrugation in the 
adsorbate-graphite interaction in Steele's expansion 
(U\). Note the sharp features of (Ulj) in both tran- 
sitions, and the sharp decrease in magnitude of (U\) 
at the commensurate- incomcnsuratc/hcrringbonc- 
ncmatic transition. 



reduction" can be significant. (For more com- 
ments on the stiffling of tilting by the substrate 
potential see Sec. IIIF.) 

B. Energy and Spreading Pressure 

In addition to structural indicators shown be- 
fore (the various order parameters), thermody- 
namic quantities are helpful in tracking phase 
transitions. Figure 3 shows the temperature de- 
pendence of the average Lennard-Jones interac- 
tion energy (Ulj) [Eq. (2)] and of the average 
corrugation potential energy (U\) [Eqs. (4) and 
(5)]. The former is useful in delineating the sys- 
tem's structural behavior irrespective of the sub- 
strate's potential, whereas the latter is a good 
indicator of atomic order of the substrate and 
therefore complements the molecular informa- 
tion provided by OP comm (see Sec. Ill A). 
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The average interaction energy is very closely 
related to how close together the hexane 
molecules are [in fact, (Ulj) oc JuLj(r)g(r) rdr, 
where g(r) is the atomic pair correlation function 
(see Sec. HIE)]. The large increase in the magni- 
tude at T\ = 138 ± 2 K is a clear indication that 
molecules in the incomensurate nematic phase 
have managed to become closer to the minimum 
of the Lennard-Jones interaction. At the melting 
transition (T% = 176 ± 3 K) molecular distances 
again become less optimal and this average en- 
ergy rises. Both of these features are evident 
in the plots of g C om(r), see Fig. 11. The average 
corrugation potential energy (U\) drammatically 
decreases in magnitude at Ti, a clear indication 
of the lack of registration of the adsorbate film 
and the substrate. 

Other energy ensemble averages have inter- 
esting features, but considerably less pronounced 
than those for (Ulj) and (U\) (typically they 
show changes in their derivatives with respect to 
T, but no sudden jumps). One conclusion that 
can be drawn from the behavior of these energy 
averages is that the herringbone/commensurate 
to nematic/incomensurate transition is driven 
by a competition between the intermolecular 
forces with the substrate corrugation (see Sec. 
IIIF below for more discussion on this). 

The specific heat at constant area ca (not 
shown) exhibits some scatter, but does present 
two significant sharp peaks at ca. 138 K and 176 
K and there is some indication that there may be 
a latent heat associated with both transitions. 
However, the energy distributions (not shown) 
near the transitions do not show any striking 
bimodality which would seem to suggest that the 
transition is not first-order. As it is usual with 
finite-size simulations, the determination of the 
order of the transitions remains elusive. 

Another quantity of interest is the spread- 
ing pressure (surface "tension") <3?. The virial 
expansion 
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FIG. 4: Relative spreading pressure Z of hexane ad- 
sorbed on graphite [see Eq. (14)]. Note that Z ~ 
for T < T\ as expected for a solid (small spreading 
pressure), and that Z approaches the ideal gas value 
(Z = 1) for T > Ti. The value of Z for the inter- 
mediate temperature region (Ti < T < T2) indicates 
that the nematic is a fluid phase. 



relates it to the interactions between molecules 
or equivalently to the center of mass pair cor- 
relation function g com (r) (see Sec. HIE). Here 
N/A is the areal density of the molecules. The 
"relative spreading pressure" Z is thus a good 
indicator of how close the equation of state is to 
that of an ideal gas (Z ~ 1) or that of a solid 
(Z ~ 0). Figure 4 shows Z for the hexane sys- 
tem. It is evident from our results that below 
T\ the system has a very small spreading pres- 
sure as expected for a solid phase, whereas both 
the nematic (T\ < T < T2) and isotropic phases 
(T > T2) are very very close to an ideal gas. 
These results are indicative of the liquid nature 
of the intermediate nematic phase. 

C. Distributions 

In conjunction to the order parameters men- 
tioned above, various distributions can also help 
illustrate the properties of the system, especially 
near the phase transitions, and in particular re- 
garding the behavior of various molecular de- 
grees of freedom, both positional, orientational 
and internal. In Figs. 5, 7, 6, and 8 we present 
the atomic height P(z), azimuthal angle P(4>), 
internal roll-angle P(ijj) and tilt angle P(0) dis- 
tributions for various temperatures. For P(z) 
we considered < z < 12 A sampling every 
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FIG. 5: Atomic height distributions P(z) for the her- 
ringbone solid at T = 130 K (solid black line); ne- 
matic at T = 160 K (red dotted line); and liquid at 
T = 200 K (green dashed line). The insert shows 
all results obtained from T = 100 to 200 K back to 
front respectively. The thin black lines represent the 
herringbone (T < 138 K), the thick black lines the 
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FIG. 6: Molecular roll angle distributions P(ip) ■ 
Please refer to Figure 5 for format information. Note 
how sharp the change in P(ip) is at the HB to ncmatic 
transition. 



Az = 0.0056 A. For P(fa) and P(ip) and P{6) we 
considered angles the [0°,180°] domain {[0°,90°] 
for P(9)}, sampling every 1°. The azimuthal an- 
gle for a molecule, fa is defined to be the angle 
formed between the direction of the smallest mo- 
ment of inertia projected onto the xy-plane and 
the x axis. The tilt angle 6 is defined to be the 
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FIG. 7: Molecular azimuthal angle distributions 
P{4>)- Please refer to Figure 5 for format informa- 
tion. Note how rapidly the position and shape of the 
peaks change at the HB to ncmatic transition, and 
how fast the peaks disappear in the nematic to fluid 
transition. 
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FIG. 8: Molecular tilt angle distributions P{6). 
Please refer to Figure 5 for format information. 

angle formed between the direction of the small- 
est moment of inertia and the z axis (9 = 90° 
corresponds to a molecule laying flat on the sub- 
strate). The internal roll angle tp is defined for a 
particular bond involving any three consecutive 
atoms as 



cos 




x iv 



<*)]■ 



x ( ri . 



(15) 



The roll-angle has a value of when the plane 
three consecutive atoms in a molecule is flat 
(parallel to the graphite basal plane) and it 
equals 90° when the plane is perpendicular to 
the graphite 21 . 
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FIG. 9: Dihedral angle distributions P(<t> d ) for 100 < 
T < 200 (back to front). Please refer to Figure 5 for 
format information. 

The atomic height distribution (Fig. 5) shows 
that essentially no molecules have been pro- 
moted to the second layer in the low temper- 
ature phase, and only a very small fraction in 
the nematic or liquid states. The double peak 
stucture at low temperatures can be interpreted 
easily upon consideration of the molecular roll- 
angle distributions (Fig. 6): as many molecules 
are rolled "on their side", raising some atoms 
of the hexane chain above the main peak. At 
higher temperatures molecules are flatter against 
the surface and the double-peak in P(z) disap- 
pears. The azimuthal angle distribution (Fig. 
7) shows three distinct temperature regions cor- 
responding to the HB phase for T < 138 K 
\P{4>) peaked at 30° and 130°], a nematic for 
(138 K < T < 176 K) [P(4>) peaked at 0°, but 
note the presence of a small peak near 90° which 
corresponds to the formation of domains], and 
a nearly isotropic fluid for T > 176 K. A dis- 
cussion of the tilt angle distribution (Fig. 8) is 
presented later (Sec. IV D), as it is linked to the 
mechanisms for the phase transitions. See also 
the discussion on orientational order parameters 
in Sec. Ill A. 

D. Internal degrees of freedom 

It is important to monitor the degree to which 
the phase transitions exhibited by the system are 
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FIG. 10: End-to-end distance distributions P(d) for 
100 < T < 200 (back to front). Please refer to Figure 
5 for format information. 



correlated to internal degrees of molecular free- 
dom. The bond-angle distribution (not shown) 
is a simple gaussian corresponding exactly to the 
Maxwell-Boltzmann distribution for the poten- 
tial Ub enc i [Eq. (6)]. This is, in fact, indicative 
that the bond angles are unaffected by interac- 
tions with other molecules or substrate and irrel- 
evant to the properties of the phases and phase 
transitions studied. Figures 9 and 10 show the 
distributions of the dihedral angles, P(4>d), as 
well as end-to end distances (CH3-CH3 separa- 
tions in the same molecule), P(d), respectively. 

The information on both figures is signifi- 
cantly correlated, as the lower-length molecules 
correspond to the presence of one or more gauche 
defects [the width of the main lobe in P{d) cor- 
responds to the width of the bond-angle distri- 
bution]. We find that the number of gauche 
defects is significantly smaller than for hexane 
molecules in vacuo [equipartition based on utors , 
Eq. (7)], mostly due to the "flattening effect" of 
the substrate potential, but also due to the in- 
teractions with the neighbors to a lesser degree. 
The number of gauche defects is small until close 
to T2 and does not present a very significant 
temperature dependence. The end-to-end prob- 
ability distribution shows a slightly more pro- 
nounced change at T%. Given that the fraction 
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FIG. 11: Center of mass molecular pair correlation 
function g C om{r) (top panel) and molecular orien- 
tational pair correlation function g<j,(r) for the her- 
ringbone solid at T — 130 K (solid black line); ne- 
matic at T = 160 K (red dotted line); and liquid at 
T = 200 K (green dashed line) . The left inserts show 
the transitional region between herringbone and ne- 
matic (130 K < T < 150 K), the right inserts show 
the transitional region between nematic and liquid 
(170K<T< 180 K). 



of molecules involved (those with gauche defects) 
is so small, we find it unlikely that these internal 
degrees of freedom affect significantly the phase 
transition at T\. For the phase transition at T2 
we find that gauche defects do contribute to the 
footprint reduction as first suggested by Hansen 
et a/. 4,5 , however that contribution does not ap- 
pear to be as important as that from tilting (see 
Sec. IIIF). 

E. Correlations 

We have also obtained various correlation 
functions that permit a more complete elucida- 
tion of the characteristics of the three phases ob- 
served for hexane on graphite and to calculate 
properties that can be directly probed experi- 
mentally [i.e. the static structure factor S(q), see 



below]. 

We start by defining the center of mass molec- 
ular pair correlation function g C om( r ), which cor- 
responds to the joint probability of finding a 
molecule with center of mass a distance r away 
from another. We also define the orientational 
pair correlation function fl^(r), which is calcu- 
lated by averaging cos2(c/>j — fa) (fa is the az- 
imuthal angle of the i-th molecule) calculated 
over all molecule pairs (i,j) whose centers of 
mass have a separation r. (In practice correla- 
tion functions are calculated on a discrete sam- 
pling of the distance between the centers of mass, 
in our case we used Ar = 0.083 A, also note that 
we use two-dimensional normalization.) 

Figure 11 presents g com (r) and g^r) for 
various relevant temperatures. For 130 K, 
both curves confirm that the HB phase is solid 
[9com(r) has a sharp peak structure and vanishes 
between peaks], and the sharp alternating signs 
of g^>(r) is indicative of the HB phase too. At 
160 K, however, the correlation functions reveal 
that the mesophase is, indeed, liquid, but with a 
very strong orientational order corresponding to 
a nematic phase. The evidence is quite strong to- 
wards labeling the mesophase as a nematic liquid 
crystal. At 200 K, g CO m(r) is rather smooth and 
g^(r) rapidly decays to 0, a clear indication of 
a regular liquid behavior. This behavior is con- 
sistent with the conclusions obtained from the 
analysis of the spreading pressure in Sec. IIIB 
(see Fig. 4). 

A second important correlation function cor- 
responds to the atom-atom pair correlation func- 
tion g(r) which is defined as the joint proba- 
bility of finding one atom at a distance r from 
another (again, distances were discretized with 
Ar = 0.083 A, and we used two-dimensional nor- 
malization). Figure 12 shows the atomic pair 
correlation function for three characteristic tem- 
peratures corresponding to a HB (130 K), a 
nematic (160 K) and a liquid (200 K). In all 
cases, we observe very sharp peaks correspond- 
ing to intramolecular pairs: r = 1.54, 2.58, 3.96, 
5.16 and 6.50 A which correspond exactly to the 
atomic distances in a hexane molecule in its low- 
est configurational state. Broader (and weaker) 
peaks (visible only in the lower panel) at small 
distances correspond to intramolecular pairs af- 
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FIG. 12: Atomic pair correlation function g(r) for 
the Herringbone solid (top panel, T = 130 K), nc- 
matic (middle panel, T = 160 K), and liquid (bottom 
panel, T = 200 K). The sharp peaks numbered 1-5 
correspond to intramolecular pairs at distances cor- 
responding to a hcxane molecule in its lowest config- 
urational state: r = 1.54, 2.58, 3.96, 5.16 and 6.50 A. 
Broader (and weaker) peaks (visible only in the lower 
panel) at small distances correspond to intramolecu- 
lar pairs affected by a gauche defect: r = 3.11, 4.55 
and 6.05 A. The broad and relatively smooth behav- 
ior at larger distances corresponds to pairs belonging 
to different molecules. 



fected by a gauche defect: r = 3.11, 4.55 and 
6.05 A. The broad and relatively smooth behav- 
ior at larger distances corresponds to pairs be- 
longing to different molecules. 

From the atomic pair correlation function, 
one can readily obtain the static structure factor 
S(q) which is very important because it contains 
spectral weight information in reciprocal space, 
and is therefore directly connected to experimen- 
tal diffraction data 2-6 . Since most experiments 
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FIG. 13: Static structure factor S(q) for the Herring- 
bone solid (top panel, T = 130 K), nematic (mid- 
dle panel, T = 160 K), and liquid (bottom panel, 
T = 200 K). 



are performed with powder samples, we chose to 
calculate a three-dimensionally averaged static 
structure factor: 



5(g) = 1 + 



sin qr 
qr 



[g(r) — l]r dr . (16) 



Figure 13 presents the static structure factor cal- 
culated from the atomic pair correlation func- 
tions in Fig. 12. We note a general broadening 
of S(q) with increased temperature and only a 
rather minor shift towards smaller momenta of 
the main peak. The results are in qualitative 
agreement with experimental evidence 2-6 . 
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F. Variations 

In order to better understand the nature of 
the phase transitions (in particular the less well 
known HB to nematic) we ran MD simulations 
where some of the interaction potentials were 
varied. We considered three interesting possi- 
bilities: (a) variations of the Lennard- Jones in- 
teraction parameter ecH 2 an d £ch 3 correspond- 
ing to the pseudo atoms of the hexane molecules 
[Eq. (2)], from 25% to 200% of its "standard" 
value; (b) variations of the corrugation potential 
strength in Steele's expansion [Eqs. (4) and (5)], 
from 25% to 1,000% of the commonly accepted 
strength. ; and (c) variations of the dihedral po- 
tential strength u tor s{4>d) [Eq. (7)] from its usual 
strength (Table II) to ten times higher thus ef- 
fectively eliminating gauche defects. 

We find that for the first type of variations 
[item (a) above], both T\ (the HB to nematic 
transition) and Ti (the nematic to isotropic liq- 

1/2 

uid transition) scale roughly like e CH with a 
correlation R 2 ~ 0.999. As T2 — T\ increases 
for larger ecH n the nematic has a better chance 
of "expressing itself" and the maximum value 
of OP nem increases too. This simple scaling 
indicates that the substrate-adatom interaction 
[mostly through the Eoi(zi) term — see Eq. (4), 
which holds the molecules stronger to the sur- 
face stifling tilting] is dominant in setting the 

temperature scales [since Lorentz-Bertholt rules, 

1/2 

Eq. (3), makes this interaction scale like e CH ]. 
The reason the holding potential influences the 
phase transitions directly can be attributed to 
its stiffiling effect on tilting and the formation of 
gauche defects. 

In the second type of variations [item (b) 
above], we find that T\ increases with the corru- 
gation energy, whereas T2 is almost insensitive 
to it. As T\ increases, the nematic phase be- 
comes "weaker", eventually disappearing for a 
corrugation ca. 5 times the standard value, and 
we see a direct transition between the HB and 
the isotropic fluid. These observations demon- 
strate that the corrugation is irrelevant for the 
nematic to isotropic melting (as expected, since 
this transition occurs between two incommensu- 
rate phases). Furthermore, the increase of T\ 
is an indication that, in addition to the effect of 



holding potential described above, a competition 
between the corrugation potential (which favors 
a commensurate phase) and the interactions be- 
tween molecules (which favors some other char- 
acteristic intermolecular spacing) is important. 

The third type of variation [item (c) above] is 
useful to characterize which degrees of freedom 
contribute to the footprint reduction mechanism 
that creates in-plane space that facilitates the 
phase transitions. Since the transition at T\ hap- 
pens when the number of gauche defects is mi- 
nuscule, it is evident that they are not important 
there. For the phase transition at T2 we found 
(Sec. Ill) a significant increase in tilting and 
gauche defects as the system achieves its high 
temperature phase. Our variation of ut ors (0rf), 
which essentially eliminates gauche defects, re- 
sulted in an increase of T2 by ca. 20 K. This 
clearly indicates that gauche defects are impor- 
tant for the footprint reduction as first suggested 
by et al. ' 5 , however tilting appears to be more 
important. 

IV. DISCUSSION 

A. Herringbone commensurate solid to 
nematic liquid crystal transition at 
Ti = 138 ±2K 

The low temperature (T < 138 K) phase of 
the hexane on graphite system is a commen- 
surate herringbone solid, as evidenced by the 
structure shown in Fig. 1, the values of OPherr, 
OPcomm and OP nem in Figure 2, and the two dis- 
tinct peaks in the azimuthal angle distribution 
P(<t>) at 30° and 150° in Fig. 7. The spread- 
ing pressure of this solid system is very small 
(low temperature region of Fig. 4). Since the 
the monolayer is complete the system creates in- 
plane room by the molecules spending a signifi- 
cant amount of their time rolled on their side, as 
seen in the large peak in the roll angle distribu- 
tions P(4>) in Fig. 6 as well as the double-peaked 
atomic height distribution P(z), see Fig. 5. 

As the temperature of the system is in- 
creased, the usual signatures of increasing ther- 
mal fluctuations appear, up until T\ = 138 ± 2 
K, where the system undergoes a phase change. 
Together, the structural order parameters in Fig- 
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ure 2 indicate that the system becomes incom- 
mensurate at the same time it loses its herring- 
bone orientational ordering. In fact, this is also 
confirmed by the sharp decrease in magnitude 
in the average Steele corrugation energy (U\) 
(Fig. 3, bottom panel) indicates that, even for 
atomic degrees of freedom the system is incom- 
mensurate. The center of mass pair correlation 
function g com {r) (Fig. 11, top panel) indicates 
that the new phase is fluid, as evidenced by the 
fact that g C om(f) has relatively smooth features, 
and that the molecules are slightly closer to one 
another, which also leads to an increase in the 
magnitude of the interaction energy {Ulj) (Fig. 
3, top panel). The liquid nature of this phase is 
also evidenced by the relatively large spreading 
pressure (Fig. 4). As the solid-nematic transi- 
tion proceeds, the tilt angles of the molecules 
fluctuate creating a footprint reduction that fa- 
cilitates the phase transition (yet tilting is once 
again reduced when the new phase is reached), 
as seen by the sharp peak in OPtut at ca. 138 K 
(Fig. 2, bottom panel). 

Inspection of the snapshots in Fig. 1, plus 
more quantitative considertions based on OP nern 
and P{4>) further indicates that the positional 
order lost in this phase transition occurs simul- 
taneously with a new orientational order where 
all molecules have a common director. The 
molecules experience, in fact, a relatively long 
range orientational order (see 160 K curve, Fig. 
11, bottom panel). We believe that this phase 
has all the characteristics of a two dimensional 
nematic liquid crystal. Analysis of the probabil- 
ity distributions for the bond and dihedral an- 
gles, and for the molecular end-to-end distance 
indicate that the internal degrees of freedom do 
not cooperate significantly with the HB to ne- 
matic phase transition. In fact, it seems that 
they hardly participate at all in the relevant pro- 
cesses in the temperature range studied for a 
monolayer of hexane adsorbed on graphite. It 
should be noted that number of gauche defects 
is significantly smaller than what would be pre- 
dicted in vacuo, which is due to a large degree to 
the substrate interaction, and to a lesser degree 
to the interaction between molecules which has 
an aligning effect. 



B. Nematic to isotropic liquid phase 
transition at T 2 = 176 ±3K 

As the temperature is raised, the nematic 
persists up until Ti = 176 ± 3 K where another 
orientational phase transition occurs, this time 
to an isotropic fluid. The phase transition is 
characterized by a sudden drop in the nematic 
order parameter OP nem (Fig. 2, center panel) 
and the disappearance of any directional pref- 
erence in the azimuthal angle distribution P((f>) 
(Fig. 7, see also the random orientation in the 
200 K snapshot in Fig. 1). In addition, the ori- 
entational pair correlation function g<f,{r) shows 
that angular correlations between molecules de- 
cay extremely rapidly with distance in the high 
temperature phase (see 200 K curve, Fig. 11, 
bottom panel). The Lennard- Jones interaction 
increases sharply (which is related to the in- 
crease in the typical intermolecular distances), 
and there is a strong peak at at Ti in the spe- 
cific heat ca- 

In the liquid, the tilt angle behavior shown 
in the bottom panel of Fig. 2 reveals that the 
molecular axes are much freer to fluctuate than 
in the nematic and this may cause a signifi- 
cant reduction in the molecular footprint (as 
discussed in Sec. IIIF, gauche defects also con- 
tribute to this footprint reduction, though tilting 
appears to be more important). The dihedral 
angle and end-to-end distance distributions in 
Figs. 9 and 10 show that, in the liquid, gauche 
and other dihedral defects are promoted more 
readily with increasing temperature, but, as with 
the solid-to-nematic phase transition, no new or 
abrupt changes in the internal degrees of molec- 
ular freedom accompany melting. 

C. Information related to scattering 

The atomic pair correlation functions in Fig. 
12 are a combination of intra- and intermolec- 
ular pair terms and hence are rich in infor- 
mation about how the internal degrees of free- 
dom behave as well as the bulk structure of 
the system's phase. The five extremely sharp 
peaks correspond to atom-atom distances within 
the molecules, and correspond precisely to the 
intramolecular atomic distances in the ground 
state of hexane. Careful inspection of the 
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g(r) curves at high temperatures reveals the 
emergence of smaller more diffuse peaks, which 
again precisely correspond to the intramolecular 
atomic distances, this time in presence of one 
gauche defect. When the curve g(r) starts to 
elevate off the horizontal axis, the contribution 
from atom pairs between different molecules be- 
gins and the decreasing structural order with 
increasing temperature is evident. The static 
structure factor S(q) shown in Fig. 13 helps fur- 
ther understand the behavior of g(r). We note a 
general broadening of S(q) with increased tem- 
perature and only a rather minor shift towards 
smaller momenta of the main peak. The results 
are in qualitative agreement with experimental 
evidence 2-6 . 

D. Mechanisms 

As in previous computer simulations 7-9 , our 
work reproduces general features of the real sys- 
tem as explored by experiment, and quite accu- 
rately determines the melting temperature T2 1 " 6 . 
Our findings of an intermediate nematic liquid 
crystalline phase between T\ and T2 are roughly 
in agreement with previous simulations by Pe- 
ters and Tildesley 8 . However, the existence of 
this intermediate phase has not been confirmed 
experimentally, which suggests that larger com- 
puter simulations would be of interest in deter- 
mining the presence of highly complicated do- 
main behavior, or even dynamic domains that 
smaller simulations simply cannot capture. 

A unique benefit of computer simulations is 
that they allow careful inspection of the mecha- 
nism for various phase transitions when enough 
data is present. In the work presented here, 
we see that in the solid herringbone phase the 
molecules spend much time rolled so that their 
backbones are ziz-zagged, or rough relative to 
the graphite substrate (see Fig. 6). The solid 
to nematic transition involves negligible amount 
of gauche defects (Fig. 9), instead showing pro- 
nounced molecular tilt fluctuations (Fig. 2, bot- 
tom panel) in concert with the molecules' rolling 
on their sides. The resulting nematic phase 
still has an average density of unity but the 
molecules are closer to each other (in the di- 
rection perpendicular to the molecules), as ev- 
idenced by the shift in the first peak of the cen- 



ter of mass pair distribution function gcorrSj^) 
(Fig. 11) and the sudden increase in magnitude 
of the Lennard-Jones interaction {Ulj) (Fig. 3, 
top panel). As a result, the tilt-angle distribu- 
tion (Fig. 8) takes on a fundamentally different 
shape somewhat more sharply peaked that the 
gaussian-like distributions seen in the solid. The 
importance of the molecule-molecule interaction 
in the mesophase is underscored by the fact that 
the tilt angle distributions do not considerably 
broaden in the nematic. 

At a higher temperature the system melts 
onto an isotropic phase. Although no sudden 
increase in tilt fluctuation accompanies melting, 
the tilt angle distributions begin to show ther- 
mal broadening, accounting for the increased 
slope of OPtut in the liquid (Fig. 2, bottom 
panel). A possible "signature" of the nematic 
is that the tilt angles are locked and that it 
is a "tilt-locking" mechanism that preceeds the 
mesophase. The first simulations completed on 
the system 4 ' 5 concluded that, at melting, pro- 
liferation of dihedral (gauche) defects in the 
adsorbate cause the projected presence of the 
molecules (their footprint) to decrease, thus 
creating more in-plane room for center-of-mass 
molecular fluctuations and thus allowing melting 
to ensue. In fact, freezing of the dihedral degrees 
of freedom resulted in a considerable increase in 
the system's melting temperature. Later work 
by the same authors 6 and others 7-9 suggests 
that, when the molecule-substrate interaction 
strength is decreased to more realistic values, 
guche defects become unimportant for melting. 
Our work suggests that the creation of in-plane 
free space is indeed important for all phase tran- 
sitions in the system, but that the "footprint 
reduction" happens exclusively through tilting 
(and moderate stacking) at the solid-nematic 
transition, where the internal molecular degrees 
of freedom are unimportant. The "footprint re- 
duction" for the nematic-isotropic transition is, 
indeed, caused by both gauche defects and tilt- 
ing, as both the probability of finding signifi- 
cant dihedral angles and tilting start rising sig- 
nificantly near T2 (see bottom panel of Fig. 2 
and Fig. 9). 

The simulations conducted on the variation 
of various Lennard-Jones interaction parameters 
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agree with the above assessments because when 
only the lateral graphite corrugation well depth 
is varied [type (b) as defined in Sec. IIIF], T\ 
increases slightly and T2 remains unchanged. 
However when the potential parameters for the 
pseudoatoms themselves are varied [type (a)], 
the interaction of the adatoms with the sub- 
strate also increases, including a stronger verti- 
cal force component. Thus, tilting is stifled and 
T\ and T2 are significantly affected. Since we 
already know that an increase in lateral forces 
alone cause relatively small changes in T%, such 
results underscore the importance of tilting prior 
to the nematic. Our elimination of gauche de- 
fects by varying the torsional potential [variation 
type (c)] resulted in a moderate increase of the 
transition temperature T%, from which we con- 
clude that both tilting and gauche defects are 
significant, albeit it seems that the former are 
more important. In essence, we agree qualita- 
tively with earlier work's findings that "footprint 
reduction" is important in this system, only in 
our work the prevalent mechanism is tilting. 



nematic to isotropic transition is, as expected, 
rather insensitive to the corrugation potential, 
but remains very sensitive to the main part of 
the substrate-adatom potential. 

We have observed somewhat contradictory 
evidence in terms of trying to identify the or- 
der of the phase transitions. Whereas the spe- 
cific heat curves seem to indicate that there may 
be a latent heat in both cases, the lack of bi- 
modality in the energy distributions would in- 
dicate otherwise. The inconclusiveness is not, 
however, surprising given the finite size of the 
system simulated. 

We have also performed MD simulations with 
ca. 2x and 4x the number of molecules dis- 
cussed throughout this paper. We find that the 
HB to nematic phase transition temperature T\ 
is only weakly dependent on system size (but 
the relatively poor statistics possible with the 
larger sizes make it difficult to extrapolate to the 
thermodynamic limit). The nematic to isotropic 
phase transition temperature T2 is insensitive to 
system size for the sizes considered. 



E. Final Remarks 

Our simulations indicate that the HB to ne- 
matic transition is due to a competition between 
the pro-commensurate interaction between ad- 
sorbate molecules and the corrugation substrate 
potential and the intermolecular interactions 
which favors intermolecular distances different 
from the required by the substrate. The latter 
prevail at ca. 138 K, leading to a nematic liquid 
crystalline state. If a stronger attraction to the 
substrate is used in our similations (Sec. IIIF) 
the stifling of mechanisms that reduce the molec- 
ular footprint makes the phase transition occur 
at a higher temperature. An increase in the 
corrugation potential has a similar effect. The 
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